﻿using System;
using CloudDetector.Utils;

namespace CloudDetector.Algorithms
{
    /// <summary>
    /// Otsu küszöbölést megvalósító osztály
    /// </summary>
    public class Otsu
    {
        #region Fields

        private ImageLayers Source;

        #endregion

        #region Contructor

        public Otsu(ImageLayers S)
        {
            Source = S;
        }

        #endregion

        #region Member Method

        /// <summary>
        /// Otsu köszöbölés minden rétegen egyszerre.
        /// </summary>
        /// <param name="Result">Eredmény rétegek.</param>
        /// <param name="Cloud">Felhőrétegek.</param>
        /// <returns>Küszöbértékek tömbje.</returns>
        public int[] Thresholding(out ImageLayers Result, out ImageLayers Cloud)
        {

            int LayerCount = Source.Layercount;
            int H = Source.Height;
            int W = Source.Width;
            int[,] Histogram = new int[LayerCount, 256];

            //histogram elkészítése
            for (int j = 0; j < H; ++j)
                for (int i = 0; i < W; ++i)
                    for (int k = 0; k < LayerCount; ++k)
                        ++Histogram[k, (Source.layers)[k, i, j]];


            double[] Mu = new double[LayerCount];
            double[,] P = new double[LayerCount, 256];       //normalizált hisztogram
            int PixelCount = H * W;

            //hisztogram normalizálás és mű (Mu) kiszámolása
            for (int i = 0; i < 256; ++i)
                for (int j = 0; j < LayerCount; ++j)
                {
                    //ha nincs itt a (double) akkor int osztás :S
                    P[j, i] = Histogram[j, i] / (double) PixelCount;
                    Mu[j] += i * P[j, i];
                }

            //mű_k-k és q kiszámolása a rekurzív képlettel, és szigma_négyzet_b kiszámolása
            double[,] Q = new double[LayerCount, 256];
            double[,] Mu1 = new double[LayerCount, 256];
            double[,] Mu2 = new double[LayerCount, 256];
            double[,] Szig2B = new double[LayerCount, 256];

            for (int k = 0; k < LayerCount; ++k)
                Q[k, 0] = P[k, 0];
            double[] Max = new double[LayerCount];
            int[] Thresholds = new int[LayerCount];                      //köszüb értéke
            for (int i = 0; i < 255; ++i)
            {
                for (int k = 0; k < LayerCount; ++k)
                {

                    Q[k, i + 1] = Q[k, i] + P[k, i + 1];
                    Mu1[k, i + 1] = (0 != Q[k, i + 1]) ? (Q[k, i] * Mu1[k, i] + (i + 1) * P[k, i + 1]) / Q[k, i + 1] : 0;
                    Mu2[k, i + 1] = (0 != (1 - Q[k, i + 1])) ? (Mu[k] - Q[k, i + 1] * Mu1[k, i + 1]) / (1 - Q[k, i + 1]) : 0;
                    Szig2B[k, i + 1] = Q[k, i + 1] * (1 - Q[k, i + 1]) * Math.Pow(Mu1[k, i + 1] - Mu2[k, i + 1], 2);
                    if (Max[k] < Szig2B[k, i + 1])
                    {
                        Thresholds[k] = i + 1;
                        Max[k] = Szig2B[k, i + 1];
                    }
                }
            }

            //küszöbölt "kép" és felhő generálása
            Result = new ImageLayers(LayerCount, W, H);
            Cloud = new ImageLayers(1, W, H);

            for (int j = 0; j < H; ++j)
                for (int i = 0; i < W; ++i)
                {
                    for (int k = 0; k < LayerCount; ++k)
                    {
                        if ((Source.layers)[k, i, j] > Thresholds[k])
                            (Result.layers)[k, i, j] = 255;
                        else
                            (Result.layers)[k, i, j] = 0;
                    }
                    //a képpont egy felhőhöz tartozik
                    //itt végrehajtani sokkal gyorsabb mint külön fv-ből újra végig menni a pixeleken
                    if ((Result.layers)[0, i, j] == 255 &&     //br
                        (Result.layers)[1, i, j] != 255 &&     //wn     vagyis 0
                        (Result.layers)[2, i, j] == 255 )       //hz
                    {
                        (Cloud.layers)[0, i, j] = 255;
                    }
                    else
                    {
                        (Cloud.layers)[0, i, j] = 0;
                    }
                }

            return Thresholds;
        }

        #endregion
    }
}
